clear all;
%% Code explanation
% 170305 DZ 
% 

%% Define the resistance to temperature conversion equation
T_RuOx=@(R)1./exp(-1.001200439904+1.175279377873*log(R-2.2)+0.021010648114*(log(R-2.2)).^2-0.024007137972*(log(R-2.2)).^3+0.002054398372*(log(R-2.2)).^4);
%
%% Set the parameter

L_W1=121.9/105.5; % 3.3Sn/12PbTe #766
L_W2=85.2/105.5; % 3Sn/12PbTe #849

%% Read the data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


root=('D:\Documents2\Rawdata\Sn\171219 MPI mK\20171211 data 3Sn 12PbTe inplane\'); %

filenames={'20171214_3samples_baseT';'20171215_3samples_50uW';'20171215_3samples_100uW';...
    '20171215_3samples_180uW';'20171215_3samples_250uW';'20171215_3samples_350uW';...
    '20171216_3samples_200uW_noroots';'20171216_3samples_350uW_noroots';'20171216_3samples_450uW_noroots';...
    '20171216_3samples_600uW_noroots';'20171216_3samples_750uW_noroots';...
    '20171216_3samples_1000uW_noroots';'20171216_3samples_1300uW_noroots';'20171216_3samples_1600uW_noroots';...
    '20171220_3samples_2100uW_noroots';'20171220_3samples_2700uW_noroots';'20171221_3samples_3500uW_noroots'};
theta0=[7.75 7.4 7]';


NT=length(filenames);


colorset=jet(NT);
figure

Tmc=zeros(NT+1,2);
Tmc_err=zeros(NT+1,2);
Hc2=zeros(NT+1,2);

%%% 50%Rn
ratio=0.5;
Tmax1=0.938;
Tmax2=1.0895;




Tmc(NT+1,1)=Tmax1;
Tmc(NT+1,2)=Tmax2;

%% get the normal state resistance from high field data at the highest temperature
kk=NT;

    file1=[root,filenames{kk},'_inplane_measdown.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%f %f %f %f %f %f %f','commentStyle','#');
        data1=data1';
      
        fclose(fid);
    
      theta=round(data1{6}*100)/100;
      
      
    
          index=find(theta==theta0(2));
          
          B=data1{1}(index);
          R1=(data1{2}(index)/L_W1+0.009)/5*1000;
         
          indexB=find(B>4);
          p1=polyfit(B(indexB),R1(indexB),1);
           subplot(2,2,1),plot(B,ratio*(p1(1)*B+p1(2)),'-k');
          hold on
      
   
          index=find(theta==theta0(3));
          
          B=data1{1}(index);
          R2=(data1{3}(index)/L_W2+0.02)/5*1000;
       
        indexB=find(B>4);
        p2=polyfit(B(indexB),R2(indexB),1);
        subplot(2,2,2),plot(B,ratio*(p2(1)*B+p2(2)),'-k');
          hold on
          
      

%% plot the data and get Hc2
for kk=1:NT
   
    file1=[root,filenames{kk},'_inplane_measdown.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%f %f %f %f %f %f %f','commentStyle','#');
        data1=data1';
        fclose(fid);
    
      theta=round(data1{6}*100)/100;
      
     
          index=find(theta==theta0(2));
          B=data1{1}(index);
          R1=(data1{2}(index)/L_W1+0.009)/5*1000;
          T=T_RuOx(data1{5}(index)*100);
          
          subplot(2,2,1),plot(B,R1,'Color',colorset(kk,:)','LineWidth',1);
          hold on
          Tmc(kk,1)=mean(T);
          Tmc_err(kk,1)=3*std(T);
          
          index1=find((R1-ratio*(p1(1)*B+p1(2)))<=0);
        if isempty(index1)==0
          Hc2(kk,1)=B(index1(1));
        end
  
 
          index=find(theta==theta0(3));     
          B=data1{1}(index);
          R2=(data1{3}(index)/L_W2+0.02)/5*1000; 
          T=T_RuOx(data1{5}(index)*100);
          subplot(2,2,2),plot(B,R2,'Color',colorset(kk,:)','LineWidth',1);
          hold on
         
          Tmc(kk,2)=mean(T);
          Tmc_err(kk,2)=3*std(T);
          
           index2=find((R2-ratio*(p2(1)*B+p2(2)))<=0);
        if isempty(index2)==0
          Hc2(kk,2)=B(index2(1));
        end
          
       title={'B(T)','R(Ohm)','T(K)','averaged T(K)'};
        AA=[B R2 T Tmc(kk,2)*ones(length(B),1)];
       
         xlswrite('Fig2E.xlsx',AA,kk);
       xlswrite('Fig2E.xlsx',title,kk);
      
%     
end

figure
plot(Tmc(:,1),Hc2(:,1),'dr',Tmc(:,2),Hc2(:,2),'sb');
hold on



%% End




